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1  Introduction 

Stability  of  shear  flows  at  high  Reynolds  number  has  been  one  of  the  cornerstones  of  hy¬ 
drodynamic  stability.  A  canonical  problem  which  contains  the  basic  features  is  that  of  a 
shear  layer  -  two  streams  of  different  velocity  flowing  past  each  other  with  different  density. 
Helmholtz[18]  and  Kelvin[21]  were  probably  the  earliest  to  consider  evolution  of  distur¬ 
bances  in  a  stratified  vortex  sheet.  While  attempting  a  more  general  case  Rayleigh[27] 
approximated  the  shear  layer  by  piecewise-linear  velocity  profiles. 

A  kinematic  description  of  the  instability  for  the  unstratified  vortex  sheet  was  given  by 
Batchelor [1]  purely  based  on  vortex  dynamics  (figure  1(a)).  Imagine  perturbing  the  vortex 
sheet  by  a  sinusoidal  disturbance  so  that  the  perturbed  interface  is  located  at  77  =  sinkx. 
In  the  neighbourhood  of  the  nodes  (A  and  B),  the  positive  vorticity  induces  a  clockwise 
circulating  velocity  field.  If  drj/dx  >  0,  the  crest  and  trough  move  away  from  each  other, 
leading  to  vorticity  being  swept  off  from  nodes  like  A,  whereas  if  drjjdx  <  0,  the  crest 
and  trough  come  closer  to  each  other,  leading  to  vorticity  being  swept  into  nodes  like  B. 
Thus,  accumulation  of  vorticity  at  points  like  B  takes  place  unboundedly  in  the  linear, 
non-dissipative  scenario,  giving  exponential  growth.  In  absence  of  stratification  the  vortex 
sheet  problem  as  studied  by  Helmholtz  and  Kelvin  is  unstable  to  infinitely  large  wavenum¬ 
bers.  This  ill-conditioned  nature  of  the  system  gets  removed  in  Rayleigh’s  problem  with  the 
introduction  of  a  length  scale  (figure  1(b)).  Rayleigh’s  problem  consists  of  two  interfaces 
of  vorticity  discontinuity.  The  Kelvin-Helmholtz  (KH)  instability  of  a  shear  layer  was  the 


(a)  (b) 


Figure  1:  (a)  Vortex  sheet  (b)  Mixing  layer 
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Figure  2:  Kelvin-Helmholtz  instabilties 


earliest  demonstrations  of  how  inviscid  instabilities  can  occur  due  to  presence  of  vorticity 
extremum  in  the  flow. 

KH  instability  is  commonly  observed  in  nature.  In  atmospheric  context  clouds  often  help 
to  visualise  them.  In  figure  2  KH  instabilities  are  seen  in  a  laboratory  experiment  and  also 
in  rivers. 

Though  KH  instability  considers  the  two  layers  to  be  of  two  different  density,  for  stable 
stratification  the  role  of  gravity  is  purely  stabilising.  Taylor[32]  and  Goldstein[17]  simulta¬ 
neously  showed  how  stable  stratification  can  have  a  destabilising  role  when  it  interacts  with 
shear.  Considering  various  multi-layer  velocity /density  profiles  Taylor  studied  instabilities 
arising  out  of  interaction  of  vorticity  and/or  gravity  waves.  We  will  revisit  one  of  Taylor’s 
problems  later  in  the  analysis. 

Subsequently  to  model  a  stratified  shear  layer  relevant  for  geophysical  problems,  Holmboe[20] 
considered  Rayleigh’s  piecwise-linear  approximation  with  an  embedded  density  interface  at 
the  shear  layer  centre  (figure  3).  For  weak  stratification  a  KH  like  behaviour  was  observed, 
instabilities  stationary  with  respect  to  the  flow.  On  increasing  stratification  a  new  class  of 
instabilities  emerged  -  a  pair  of  counter  propagating  waves  known  thereafter  as  Holmboe 
waves.  Unlike  KH  instability  Holmboe  waves  are  not  easily  visualised  in  nature  due  to  the 
nonstationary  behaviour  of  the  waves.  They  have  been  experimentally  realised  and  consid¬ 
ered  to  be  present  in  exchange  flows.  Often  a  distinction  is  made  between  KH  and  Holmboe 
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Figure  3:  Stratified  shear  layer 


waves  based  on  the  stationarity  nature  as  one  moves  with  the  speed  of  the  shear  layer  centre. 
Such  a  classification  though  accurate  in  the  symmetrical  case  becomes  ambiguous  in  the 
asymmetrical  case.  Recently [8]  it  was  shown  that  a  classification  based  on  wave  interactions 
-  vorticity-vorticity  or  vorticity-gravity  -  provides  a  consistent  explanation. 

Approaches  towards  understanding  evolution  of  disturbances  in  stratified  shear  layers  can 
be  mostly  classified  in  two  categories  -  Linear  stability  calculations  and  direct  numerical 
simulations  (DNS).  While  doing  linear  stability  studies  if  one  restricts  to  piecewise  linear 
profiles  then  dispersion  relations,  eigenfunctions  can  be  analytically  computed.  To  model 
realistic  profiles  if  one  considers  the  smoother  counterparts  then  solving  the  formidable 
Taylor- Goldstein  equation  (TGE)  is  the  only  recourse.  For  a  complete  nonlinear  descrip¬ 
tion  beyond  linear  instability,  DNS  holds  the  key.  DNS  enables  a  complete  exploration 
of  wide  range  of  length  and  time  scales  in  flows  transitioning  towards  a  turbulent  state. 
Having  said  that  the  computational  cost  and  effort  involved  in  a  complete  DNS  severely 
limits  an  exhaustive  study  in  nonlinear  dynamics  of  stratified  shear  layers,  many  questions 
go  unanswered.  How  does  increasing  stratification  alter  the  nonlinear  states  in  stratified 
shear  layers?  Does  the  flow  attain  a  nonlinear  stationary  state?  Gravity  waves  interact  in 
presence  of  shear  to  give  rise  to  Taylor- Caulfield  instability,  how  does  the  nonlinear  state 
evolve?  To  answer  many  such  questions  full  nonlinear  simulations  are  needed.  Is  there  a 
reduced  description  of  the  governing  equations  which  captures  the  essence  of  full  numeri¬ 
cal  simulations  while  giving  a  stronger  analytical  handle  on  the  physics?  The  goal  of  the 
present  work  is  to  outline  a  reductive  perturbation  theory  for  stratified  shear  layers.  The 
methodology  will  rely  on  a  previously  used  technique,  the  Aorticity  defect’  method. 

2  ‘Defect’  theory 

Inviscid  shear  flow  instabilities  occur  in  presence  of  vorticity  maximum,  e.g.  -  KH.  For 
KH  like  flows  one  can  immediately  notice  that  the  base-flow  quantities  vary  across  a  very 
small  region.  This  is  also  true  for  stratified  shear  layers  where  besides  vorticity,  variations 
of  density  occur  over  a  small  region.  To  obtain  a  reduced  description  we  will  try  to  exploit 
the  smallness  of  the  shear-layer  thickness.  The  perturbative  technique  to  be  adopted  will 
in  spirit  be  akin  to  boundary  layer  theory.  We  will  have  an  ‘outer’  region  where  the  pertur- 
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Figure  4:  Velocity  and  density  profiles  with  ‘defects’ 


bations  will  be  assumed  to  obey  linear,  steady  and  inviscid  dynamics  and  an  ‘inner’  region 
where  nonlinearity,  dissipation  and  unsteadiness  will  be  important.  A  key  assumption  of  the 
analysis  will  be  considering  the  base-flow  quantities  (velocity,  density)  to  be  monotonically 
varying  quantity  with  sharp  variations  or  ‘defects’  of  width  e  at  the  shear  layer  centre. 

While  studying  wake  behind  a  airfoil.  Gill [16]  modelled  the  base-state  to  be  a  Couette  flow 
with  slight  distortions.  This  severely  simplified  the  linear  stability  calculations  and  provided 
integral  dispersion  relationships.  Next  Lerner  and  Knobloch[24]  performed  long-wavelength 
stability  studies  on  a  distorted  Couette  flow.  The  slight  modification  to  Couette  flow  was 
assumed  to  arise  due  to  finite  amplitude  perturbations.  But  it  wasn’t  until  the  works  of 
Balmforth  and  co-workers  that  the  ‘defect’  model  -  a  Couette  flow  with  slight  distortions 
-  got  a  strong  theoretical  foundations  and  its  immense  applications  appreciated.  Using 
an  analysis  similar  to  critical  layer  theory  del-Castillo  Negrete  et  al.[15]  (Balmforth  et  al. 
[3])  showed  that  for  inviscid  shear  flows  one  can  reduce  2D  Euler  equations  to  a  nonlin¬ 
ear  integro-differential  equation  for  the  ‘defect’  vorticity.  Exploiting  the  ‘defect’  equations 
analogy  with  the  Valsov-Poisson  equation^  explicit  dispersion  relations  were  obtained  for 
various  distortion  imposed  on  the  Couette  flow.  Balmforth  and  Young [5]  went  ahead  and 
highlighted  evaluation  of  nonlinear  stationary  states  from  the  ‘defect’  equations.  The  lin¬ 
ear  viscous  problem  was  addressed  by  Balmforth  [6]  and  the  intriguing  connection  between 
viscous  and  inviscid  eigensolutions  was  discussed  in  context  to  various  ‘defect’  profiles.  The 
applicability  of  the  approach  for  visco-elastic  flows  has  also  been  studied  [4].  In  the  present 
work  we  will  try  to  obtain  better  understanding  of  stratified  shear  layers  aided  by  the  ‘de¬ 
fect’  theory. 

Consider  a  Bousinessq  fluid  of  density, 

p{x,y,t)  ^  p^  + 5p{x,y,t)  (1) 


The  evolution  equations,  vorticity-buoyancy  formulation,  in  the  Bousinessq  approximation 
can  be  written  in  the  following  non-dimensional  form. 


du  a(^,cj)  dB  1  o 

dt  d{x,  y)  dx  Re 

dt  d{x,  y)  ReSc 


(2) 

(3) 


^Vlasov-Poisson  equation  describes  time  evolution  of  the  distribution  function  of  plasma  consisting  of 
collision-less  charged  particles  in  the  zero  magnetic  field  limit. 
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The  horizontal  and  vertical  coordinates  have  been  non-dimensionalised  by  a  characteristic 
length-scale  Lq  and  time  by  Uq/Lq^  Uq  being  a  velocity  scale.  Stream- function  is  non- 
dimensionalised  by  UqLq.  (jJ  —  and  B  —  —g6p/ pruLo/U^  are  the  non-dimensional 

vorticity  and  buoyancy  respectively.  Reynolds  number,  Re  =  U^L^/v  {p  is  the  kinematic 
viscosity),  denotes  a  ratio  of  inertial  to  viscous  forces  and  Schmidt  number,  Sc  —  vjK  {k.  is 
the  mass  diffusivity),  provides  a  ratio  of  momentum  to  mass  diffusivity.  We  will  focus  on 
S^c  ^  oo,  the  relevant  regime  for  stratified  shear  layers^. 

As  evident  in  figure  4,  velocity  is  assumed  to  be  a  localised  ‘defect’  imposed  on  a  Couette 
flow  and  similarly  buoyancy  remains  largely  uniform  but  for  a  sharp  variation  at  the  ‘defect’ 


level.  Thus  if  the  variations  occur  over  a  length-scale  e,  we  can  write, 

'i!  = +  €^'ip{x,y,t)  (4) 

B  =  bm  +  e^b  (5) 

Hence  the  vorticity-buoyancy  equation  can  be  reduced  to, 

(6) 

bt  +  ybx  +  e^{ipxby  -  '>pybx)  =  {ReSc)~^V‘^b  (7) 


The  above  equations  can  be  solved  using  matched  asymptotic  theory  by  splitting  the  flow 
into  an  exterior  region  where  base  flow  quantities  remain  uniform  or  vary  monotonically 
and  an  interior  region  corresponding  to  the  ‘defect’. 

2.1  Outer  solution 

The  dynamics  in  the  outer  region  is  assumed  to  be  steady,  inviscid  and  linear  in  the  leading 
order.  Another  important  assumption  made  regarding  the  outer  solution  is  that  of  weak 
stratification.  It  is  assumed  that  the  outer  solution  for  in  the  leading  order,  is  unin¬ 
fluenced  by  buoyancy.  As  would  become  evident  that  in  doing  so  we  deal  with  a  weaker 
singularity  in  the  outer  problem.  An  assumption  of  these  form  becomes  unphysical  for 
strongly  stratified  flows,  scenarios  beyond  the  scope  of  stratified  ‘defects’. 

A  regular  perturbation  expansion  of  the  solution  in  the  outer  region  provides, 

V’  =  V’o  +  eV’i  +  •  •  •  (8) 

b  —  ebi  +  . . .  (9) 

Thus  the  leading  order  dynamics  is  governed  by  the  following  pair  of  equations. 

=  0,  ybia;  =  0  (10) 

^  V^V’o  = -2^(x,f)(5(y)  (11) 

Thus  the  outer  solution  behaves  like  an  irrotational  flow  being  forced  by  a  vortex  sheet,  a 
6  function  singularity,  at  the  ‘defect’  level  -  y  —  0  (buoyancy  also  gets  forced  by  a  density 

^In  the  context  of  passive  scalar  turbulence  Sc  ^  oo  is  known  as  the  Batchelor  regime.  Such  a  regime 
allows  for  rapid  variations  in  passive  scalar  field  due  to  stirring  by  a  smooth  velocity  held. 


359 


sheet).  This  necessitates  a  jump  in  tangential  velocity  (while  maintaining  continuity  of 
normal  velocity), 

[’/’0.lS-  =  -2^(^.‘).  W;3-=0  (12) 

As  will  be  seen  this  jump  in  the  outer  solution  will  be  smoothened  out  over  a  region  e  by  the 

inner  solution.  If  we  had  allowed  for  leading  order  buoyancy  effects  in  the  outer  vorticity 

equation,  the  ffow  would  have  been  forced  by  a  vortex  dipole  instead  of  a  vortex  sheet.  The 
consequence  of  these  stronger  singularity,  a  5'  singularity,  would  have  then  resulted  in  a 
jump  in  normal  velocity.  Weak  stratification  allows  us  to  overlook  such  a  possibility  in  the 
present  analysis. 

Introducing  Fourier  transform  of  the  streamfunction  as, 

/oo 

'tp{x,y,t)e~'''"^dx  (13) 

-OO 

equation  11  can  be  written  as. 


Thus  one  obtains. 


i’oyy  -  k'^i’Q  =  -2A{k,t)S{y) 


Q{k)A{k,t) 


(14) 


(15) 


where,  T  ^  denotes  the  inverse  Fourier  transform  and  the  Green’s  function,  G{k)^  can  be 
computed  for  both  unbounded  and  bounded  domains  (|^|  ^  1)  as. 


G{k) 


|^sech|A;|sinh|A;|(l  -  |^|), 


bounded 

unbounded 


(16) 

(17) 


Henceforth  all  analysis  will  be  presented  for  the  unbounded  case  but  the  results  have  been 
extended  to  the  bounded  case  too.  The  ‘defect’  analysis  is  not  solely  restricted  to  the  outer 
ffow  being  Couette  and  can  be  extended  to  more  general  ffows  by  the  method  of  constructing 
Green’s  function  for  ffow  profiles  with  curvature [2]. 

Next  we  will  regularise  the  discontinuity  in  outer  velocity  field  by  evaluating  the  inner 
solution. 


2.2  Inner  solution 

Similar  to  a  traditional  boundary-layer  like  analysis,  the  inner  solution  must  contain  at  least 
one  of  the  physics  that  have  been  overlooked  in  the  outer  region  -  viscosity,  nonlinearity 
and  unsteadiness.  Working  with  our  small  parameter,  the  defect  thickness  e,  we  define  a 
stretched  vertical  coordinate,  r]  —  yje  and  a  slow  time-scale  r  =  et.  It  needs  to  be  mentioned 
that  the  method  of  analysis  bears  similarities  with  standard  critical  layer  theories  ([33], [31]). 
Beyond  the  initial  transients  the  O(e^)  nonlinear  term  acting  on  a  critical  layer  of  width 
0(e)  becomes  important  after  time  0(l/e).  As  would  become  clear  that  the  ‘defect’  theory 
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is  considerably  simpler  than  critical  layer  analysis.  We  pose  the  following  expansion  in  the 
inner  variables, 


-ip  =  'ipo{x,0,t)  +  e(t)i{x,ri,T)  +  O(e^), 

(18) 

'4’y  —  (plri  +  ^(p2ri  +  0{e^), 

(19) 

V  Ip  =  €  (plriri  -|-  '^Oxx  “t"  4’2riri 

(20) 

b  =  B{x,r],T)  +  0(e) 

(21) 

From  above  equations  it  is  evident  that  leading  order  inner  solution,  '0o(^5  0,  t),  matches  the 
corresponding  outer  solution  at  ^  =  0.  As  mentioned  before  the  discontinuity  in  tangential 
velocity  in  the  outer  solution  will  be  smoothed  by  the  inner  solution, 

/oo 

Zdrj  (22) 

-oo 

where  Z  —  is  the  leading  order  defect  vorticity.  The  velocity  jump  equals  the  vorticity 
integrated  across  the  defect. 

Substituting  the  expansion  18  in  equations  6  and  7  and  collecting  the  leading  order  terms 
one  obtains  the  ‘defect’  vorticity-buoyancy  equations, 

Z^  +  =  Bx  +  (23) 

Br  +  r]Bx  +  ^xBrj  =  (24) 

(Z^B)  0,  as  \r]\  oo 

I  roo  -I  roo 

*  =  -mL  ^  = 

where,  'H{.)  denotes  the  Hilbert  transform  of  a  variable  defined  in  terms  of  the  following 
Cauchy  principal  value  integral. 


<2«) 

and  A  =  {e^Re)~^  is  representative  of  the  ratio  of  nonlinearity  and  viscosity.  In  critical  layer 
theory  there  exists  a  similar  quantity,  Haberman  parameter,  which  denotes  competition  of 
nonlinearity  and  viscosity. 

Equations  23-25  are  the  ‘defect’  vorticity-buoyancy  equations  relevant  for  stratified  shear 
layer  studies.  During  the  course  of  reduction  we  have  retained  the  necessary  physics  due 
to  unsteadiness,  nonlinearity  and  dissipation.  Though  the  obtained  nonlinear  integro- 
differential  equations  appear  more  daunting  than  the  governing  equations  we  started  off 
with  it  would  become  clear  that  the  ‘defect’  equations  are  more  amenable  to  analytical 
and  numerical  treatments.  Similar  equations  have  been  derived  in  critical  layer  theories  for 
stratified  shear  flows  [14]. 
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3  Conservation  laws 


The  inviscid,  nonlinear  ‘defect’  eqnations  (equation  23-24  for  A  =  0)  contains  certain  con- 
served  quantities,  some  deducible  from  evident  symmetries  in  the  problem  and  others  not 
so  obvious  ones.  Lets  us  define  a  spatial  average, 


<  .  >= 


.dxdrj 


In  absence  of  dissipation  both  average  vorticity  and  buoyancy  remains  conserved, 

<Z>r^0, 


momentum  remains  conserved, 
and  energy  also  is  conserved. 


<  TfZ  >r^  0 


(27) 

(28) 


fp'  —  $ 

<  -y- — 2  -  rjB  >r=  0  (29) 

As  always  is  the  case  with  ideal  fluid  besides  the  above  ones  there  exists  an  infinitude  of 
conserved  quantities  -  Casimirs.  For  2D  homogeneous  fluid  any  functional  of  vorticity  is  a 
conserved  quantity.  Enstrophy  -  a  representative  of  vortex  stretching  and  tilting  -  is  one 
such  commonly  used  candidate.  For  stratified  flows  due  to  baroclinic  generation  of  vorticity 
enstrophy  is  no  longer  conserved.  For  stratified  ‘defects’  one  instead  has  the  any  functional 
of  buoyancy  and  product  of  vorticity  with  any  buoyancy  functional  as  the  Casimirs, 

<  G{B)  >r=  0,  <  ZF{B)  >r=  0  (30) 

The  conserved  quantities  thus  obtained  help  in  obtaining  the  stationary  states  for  the  system 
and  also  comment  on  its  nonlinear  stability.  More  importantly  it  also  helps  as  a  check  in 
numerical  computations  to  ensure  that  the  numerical  method  indeed  conserve  the  necessary 
invariants. 


4  Inviscid  linear  stability 


Earlier  it  was  mentioned  that  previous  studies  on  flow  with  slight  distortions  exploited  the 
simplicity  of  the  model  in  deriving  explicit  integral  dispersion  relations.  Conventionally 
inviscid  linear  stability  demands  solution  of  Rayleigh  equation  or  in  the  case  of  stratified 
shear  flow  Taylor  Goldstein  equation  (TGE), 


'F  _  kA  1  = 

J  ^  {U  -c)  {U-  c)2 


(31) 


where  ^  is  the  disturbance  streamfunction,  U  the  background  base-flow,  N  —  —gj pj^dpj dz 
the  BruntVaisala  (buoyancy  time-scale),  k  is  the  disturbance  wave- number  and  c  is  the 
wave-speed.  The  relative  strength  of  stratification  vis-a-vis  inertia  is  denoted  by  Richardson 
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number,  Ri  —  (A^Lo/?7o)^*  Solution  of  the  above  boundary  value  problem,  even  numeri¬ 
cally,  is  no  easy  task  and  analytical  solutions  can  be  obtained  only  for  broken-line  profiles. 
The  ‘defect’  model  circumvents  this  difficulty  with  ease. 

Consider  the  ‘defect’  vorticity  and  buoyancy  to  be  decomposed  of  a  base-state  contribution 
and  a  perturbation. 


Z{x,'q,T)^  F{'q)  +  C{x,rj,T)  (32) 

B{x,ri,T)^j  N'^{y)dy  + /3{x,ri,T)  (33) 

Substituting  the  above  expressions  in  equations  23  and  24  and  considering  the  linearised 
equations  one  has. 


Cr  T  VCx  T  —  ^X) 

(34) 

(3t  T  V^x  T  '^x^  —  0 

(35) 

- 

2f  =  —k~  /  C^dr] 

J  —oo 

(36) 

On  assuming  a  normal-mode  form,  ({x^ri^r)  —  the  eigenfunctions  can  be 

expressed  as. 


fN^ 

{y-c) 

rFr, 

(y-c) 


fN^ 
{y  —  c)2 


(37) 

(38) 


and  the  dispersion  relation, 


V{c,k;F,N)  =  2k- 


Fy  I 

Xy-c)  {y-cf 


dry  =  0 


(39) 


This  is  the  stratified  version  of  the  integral  dispersion  relation  obtained  in  [3].  V  is  an 
analytic  function  of  c  but  for  the  real- axis  where  it  has  a  branch-cut.  The  branch-cut  gives 
rise  to  the  continuous  spectrum,  a  common  feature  of  inviscid  shear  flows.  The  presence 
of  the  continuous  spectrum  contributes  to  transient  growth  in  flows  which  otherwise  are 
modally  stable  -  Orr  mechanism  being  a  well  known  example.  Sometimes  the  dispersion 
relation  instead  of  having  zeros  in  the  principal  complex  plane,  vanishes  at  points  on  dif¬ 
ferent  Riemann  sheets.  While  doing  the  initial  value  problem  this  leads  to  observation  of  a 
collective  behaviour  of  the  continuous  spectrum  as  a  single  damped  discrete  mode  -  known 
as  Landau  pole  or  quasi-mode.  This  damping,  eponymously  known  as  Landau  damping  in 
plasma  physics,  is  a  purely  inviscid  phenomenon  arising  due  to  phase  mixing.  Quasi-modes 
have  been  probed  to  understand  the  mysterious  connection  between  the  inviscid  and  vis¬ 
cous  initial  value  problem  in  the  limit  of  vanishing  viscosity.  The  present  analysis  will  focus 
purely  on  unstable  modes  and  hence  no  longer  address  the  above  intriguing  phenomena. 
We  will  now  compute  dispersion  relations  for  various  stratified  shear  layer  configurations. 
Before  proceeding  on  to  calculations  for  both  broken  and  smooth  profiles  a  note  can  be 
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made  borrowing  from  Nyquist  theory.  For  a  closed  curve  C  in  the  complex  c  plane,  the 
winding  number, 


1 

27ri 


'D'jc) 

Vic) 


dc  =  N-P 


(40) 


where  N  is  the  number  of  zeros  and  P  gives  the  number  of  poles  of  V{c)  inside  C.  Considering 
T>  an  analytic  function  of  c  the  winding  number  gives  us  the  number  of  eigenmodes.  Earlier 
works  on  ‘defect’  theory  like  [3]  and  [6]  elucidates  how  Nyquist  theory  can  be  used  effectively 
for  the  ‘defect’  model. 


4.1  Holmboe  instability 

Holmboe  instability  occurs  when  there  exists  a  density  interface  embedded  in  a  mixing  layer. 
It  is  a  popular  model  to  represent  instabilities  arising  due  to  horizontal  shearing  of  density 
interfaces.  This  classical  instability,  in  the  symmetrical  case,  exhibits  a  transition  from 
the  stationary  Kelvin-Helmholtz  instability  to  the  propagating  Holmboe  waves.  Though 
Holmboe  waves  are  inherently  difficult  to  visualise  in  nature  due  to  their  non-stationary 
behaviour  they  have  been  observed  in  laboratory  experiments  -  exchange  ffow  between  two 
basins  of  different  density  being  an  example. 


4.1.1  Broken-line  profile 

For  a  broken-line  Holmboe  profile  denoted  by  the  following  defect  vorticity  and  buoyancy 
profile, 

Fn  =  /[^(^  -  ^o)  -  Sir]  +  r]o)], 

=  dSir]) 

one  obtains  the  dispersion  relation, 

~  ^  2^  ^  ~  Vo  —  =  Cr  +  ici  (43) 

It  was  previously  mentioned  that  the  ‘defect’  theory  is  suited  to  weak  stratification  cal¬ 
culations.  Thus  the  results  obtained  are  representative  of  small  Ri  cases  in  the  complete 
problem,  d  determines  the  strength  of  stratification  and  a  qualitative  comparison  is  made 
of  contours  of  growth-rate  in  d  —  k  plane  with  the  one  obtained  from  the  solution  of  TGE  for 
a  symmetric  shear  layer.  The  basic  features  are  well  captured  by  the  present  model  (figure 
5).  One  immediate  observation  is  the  differentiation  of  stationary  (c^  =  0)  and  propagating 
instabilities  (c^  ^  0). 


(41) 

(42) 


4.1.2  Smooth  profile 


For  smooth  profiles  the  integral  dispersion  relation  can  be  evaluated  using  Cauchy  residue 
theorem.  The  simplest  profile  to  consider  would  be  a  Lorentzian. 


_^f  \  1  1 

^  7T  {t]  +  1)^  +  a‘^  (77  —  1)^  -h  ’ 

2  ^  ad  1 

TT  r]‘^  + 


(44) 

(45) 


364 


Figure  5:  Holmboe  instability,  a)  A  broken-line  profile,  b)  Contours  of  growth-rate  (ci)  in 
d  —  k  plane  (/  =  2),  c)  Growth-rate  contours  in  Ri  —  k  plane  from  a  complete  linear  stability 
calculations  (Carpenter  et  aL[8]) 


For  the  Lorentzian  profile  the  dispersion  relation  is  identical  to  equation  43  with  c  in  the 
quadratic  expression  no  longer  being  the  complex  wave-speed  but  instead  replaced  by  c  ^ 
c  +  iasgn{ci). 

A  smooth  profile  which  would  be  used  in  the  numerical  calculations  and  represents  sharper 
variations  in  flow  quantities  is  a  tanh  profile. 


-^sech^ 

2A 


y-yo'] 

A  J 

(1 

VA 


—  tanh 


(^) 


(46) 

(47) 


One  could  note  that  the  both  the  smooth  profiles  chosen  above  are  functions  which  are 
nascent  delta  functions,  i.e.  as  cr,  A  ^  0  we  obtain  the  generalised  functions  form  as  given 
in  equation  42. 

From  residue  calculus  we  can  evaluate  the  following  integrals. 


rc 

Ii{A,rio,c)  =  / 

J  —  i 

r 

22(A,r?o,c)  =  / 

J  —( 


sech^ 


2  (y-yo 

A 


iv-c) 


df]  —  27ri 


sech^ 


2  /  c  -  770 


sech^ 


2  V-Vo 


(y-c) 


A  J  .  27Ti 

-2 - dv  =  -  — 


2sech 


A 

2  f  c-yo 
A 


«  =  ^  -  c),  s  =  i(l  +  sgn(Q)) 


tanh(^^)s  +  ^4r(2)(Q,) 


where  and  are  polygamma  functions  -  =  d'^^^{\ogr{z))/dz'^^^ . 

Based  on  the  above  evaluated  integrals  the  dispersion  relation  for  the  tanh  Holmboe  profile 
can  be  computed  explicitly  as, 

A  [f  {Ii{A,vo,c)  -Xi(A,-77o,c)}  +  dX2(A,0,c)]  =0  (48) 
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Figure  6:  Holmboe  instability,  a)  A  Smooth  profile,  b)  Variation  of  (dashed)  and  q 
(continuous)  with  stratification  (d)  for  /  =  2,  A:  =  1,  A  =  1/3, 770  =  1 


Though  we  have  obtained  an  analytical  dispersion  relation  for  the  tanh  profile  the  numerical 
evaluation  of  the  integrals  is  more  efficient  than  working  with  the  polygamma  functions. 

Figure  6  shows  the  stationary  and  propagating  nature  of  the  instabilities  for  a  set  of 
parameter  values.  Cases  from  the  instability  curves  will  be  later  visited  during  numerics. 
An  important  question  that  numerics  will  try  to  answer  is  how  does  the  nonlinear  states 
change  as  one  increases  stratification  and  transitions  from  a  KH  instability  to  Holmboe 
instability. 

4.2  Taylor- Caulfield  instability 

Holmboe  instability  arose  when  gravity  waves  interacted  with  vorticity  waves.  Taylor [32] 
during  the  course  of  studying  various  broken  line  profiles  observed  that  a  stable  stratification 
composed  of  two  density  jumps  gets  destabilised  if  there  is  an  imposed  shear.  This  initially 
surprising  result  can  be  understood  as  interaction  of  gravity  waves  facilitated  by  shear. 

4.2.1  Broken-line  profile 

For  Taylor- Caulfield  instabilities  there  is  no  defect  vorticity  and  instead  one  has  variations 
of  defect  buoyancy. 


^  d[S{ri  -  r/o)  +  5{ri  +  r/o)])  (49) 

As  seen  in  the  case  of  Holmboe  instability  one  obtains  polynomial  dispersion  relations  for 
broken- line  profiles, 

-  ^2  +  0  +  ^1  -  0  =  0,  r?o  =  1,  c  =  Cr  +  ici  (50) 
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B 


(a) 


Figure  7:  Taylor- Caulfield  instability,  (a)  A  Smooth  profile,  (b)  Growth-rate  (q)  variation 
with  stratification  (d)  for  A:  =  1, 770  =  1 


4.2.2  Smooth  profile 

Similar  to  Holmboe  case  one  can  calculate  stability  characteristics  for  various  smooth  pro¬ 
files.  For  a  Lorentzian  profile, 

(77  -h  1)^  +  (77  —  1)^  -h  OL^ _  ^ 

The  dispersion  relation  is  identical  to  equation  50  with  c  in  the  quadratic  expression  no 
longer  being  the  complex  wave-speed  but  instead  replaced  by  c  ^  c  -h  7(asgn(Q). 

For  a  smooth  tanh  profile. 


TT 


7V2 


d 


sech^ 


(^) 


-h  sech^ 


(52) 


Once  again  the  analytical  form  of  the  dispersion  relation  for  can  be  written  as, 

-  ^[X2(A,r?o,c)  +X2(A,  -r/o,c)]  =  0  (53) 

Figure  7  illustrates  the  stability  characteristics  for  Taylor- Caulfield  instability  as  obtained 
from  ‘defect’  theory.  The  TC  modes  are  purely  stationary  waves. 

Though  there  exists  some  linear  stability  calculations  and  DNS  calculations  very  little  is 
known  about  the  nonlinear  states  of  Taylor- Caulfield  instabilities.  Next  using  numerical 
solution  of  ‘defect’  equations  we  will  try  to  explore  the  details. 


5  Numerical  solution 

A  significant  advantage  that  the  ‘defect’  theory  offers  over  the  complete  Bousinnessq  Navier- 
Stokes  equations  is  in  numerical  computation.  As  is  commonly  known  that  the  notorious 
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advective  nonlinearity  in  Navier-Stokes  serves  as  a  severe  bottleneck  in  full-fledged  compu¬ 
tations.  The  defect  equations  for  Sc  oo, 

Zj-  +  rjZx  +  —  Bx  +  (54) 

Br  +  r]Bx  +  ^xBr,  =  0,  (55) 

{Z^B)  0,  as  \r]\  oo 

1  POO  1  roo 

^  2d7i,  ^  =  -- zdv  m 

though  nonlinear  avoids  the  difficult  nonlinearity  5((/)i,  (/)i^^)/5(x,  77)  {(^irjrj  —  Z).  In  the 
present  form  the  ‘defect’  equation  entails  solving  a  set  of  advection-diffusion  equations. 
The  method  adopted  in  the  present  work  relies  on  exploiting  the  similarities  of  the  ‘defect’ 
equations  with  the  Vlasov-Poisson  equation  and  using  a  operator-splitting  technique. 


5.1  Operator  splitting  technique 

The  operator  splitting  technique  used  is  based  on  the  method  outlined  in  [13]  for  solution  of 
Vlasov-Poisson  equation  and  in  [7]  for  studying  critical  layers  in  2D  vortices.  This  splitting 
technique  is  also  known  as  Strang  splitting [25]. 

The  integration  over  time-step  [r,  r  -h  Sr]  is  divided  into  three  stages, 

5.1.1  Advect  in  x  for  half  time-step 

Consider  the  advection  in  x  in  the  buoyancy  equation, 

Br  +  vt3x-0  (57) 

for  a  time  interval  5r/2.  The  exact  solution  is  ;B(x,  77,  T-h5r/2)  —  B{x  —  r]ST Fourier 
interpolation  is  used  to  perform  the  x  advection  -  Bik^cf^r  -h  SrjZ)  —  ;B(A;, 77, 

For  the  vorticity  equation, 

-h  T]Zx  =  Bx  (58) 

the  solution  can  be  similarly  written  in  Fourier  space  as  - 

Z{k^  +  St/2)  =  Z{k^  rj^  ikST/2  B{k^  t  -h  Sr/2) 

The  above  calculations  in  Fourier  space  is  done  using  the  FFT/IFFT  routine  in  Matlab. 


5.1.2  Advect  in  77  for  a  complete  time-step 

Both  buoyancy  and  vorticity  obeys  identical  rj  advection  equation. 


Br  +  ^xBfj  —  0 


$ 


1 


(59) 

(60) 


Before  performing  the  advection  ^x  is  computed  by  integrating  Z  in  77,  once  again  in  Fourier 
space.  The  rj  advection  occurs  with  a  velocity  ^x  -  B[x^  ^  +  Sr)  =  B(x,  rj  —  ^xSr^  r).  Due 
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to  lack  of  periodicity  we  use  linear/cubic  spline  interpolation  in  r]  instead  of  Fourier.  The 
same  method  is  followed  for  advection  of  Z.  If  advection  happens  into  the  domain  from 
both  above  or  below  the  computational  domain  then  it  is  assigned  the  far-held  equilibrium 
value. 

5.1.3  Advect  in  x  for  another  half  time-step 

The  hrst  step  is  repeated  for  another  half  time-step,  6t/2.  One  can  combine  the  hrst  and 
third  step  and  do  two  steps  instead  -  advect  in  r]  for  6r  followed  by  advection  in  x  for 
another  Sr. 

To  account  for  diffusion  of  vorticity  (A  ^  0)  a  forward  in  time  Euler  integration  is  car¬ 
ried  out  with  diffusion  operator  in  r]  discretized  using  finite-difference  technique. 

The  computational  domain  consists  of  {x^r])  G  [0,  27r]  x  [—ijmax^Vmax]-  Vmax  is  fixed  at 
10  so  that  the  vertical  extent  is  at  least  5  times  larger  than  the  shear  layer  and  the  far-field 
is  ensured  to  remain  uninfiuenced  by  the  dynamics  in  the  shear  layer.  The  choice  of  time- 
step,  5r,  and  horizontal  discretization,  h  =  2r:  jN {N  is  the  number  of  Fourier  modes  used) 
is  determined  by  the  following  stability  condition, 

Ivmaxl'^r  ^  h  (61) 

For  the  present  work  N  has  been  chosen  to  either  256  or  512  and  the  time-step.  At  = 
1  X  10“^,  5  X  10“^.  The  vertical  discretization,  A77,  is  chosen  so  as  to  obey  the  von  Neumann 
stability  criteria, 

AAr  1 
(Ary)^  ^  2 

At]  =  0.01,0.005  is  used  in  the  calculations.  A,  an  estimate  for  momentum  diffusion,  has 
been  widely  varied  from  0.05  —  10“^.  All  numerics  reported  hereafter  are  for  A  =  10“^. 
To  instigate  instabilities  in  the  ffows  the  base-state  is  initially  seeded  with  a  disturbance 
of  the  form  coskx  of  magnitude  0.1%  of  the  base-state  maxima.  The  r]  interpolation 
has  been  done  using  linear  interpolation.  The  interpl  command  in  Matlab  has  been  used 
for  the  purpose.  For  stratified  KH  (section  5.3.1)  a  grid  resolution  of  256  x  2001  has  been 
used  with  a  time-stepping  of  5  x  10“^.  The  remaining  calculations  have  been  done  on  a 
512  X  4001  grid  with  At  =  1  x  10“^. 

5.2  Analytical  comparisons 

To  check  the  efficacy  of  the  operator  splitting  technique  a  comparison  is  made  between  the 
numerical  solution  and  analytical  solution  known  for  a  passive  scalar  advection  equation. 
The  passive  scalar  advection  equation, 

Xt  +  VXx  +  2X7?  cos  2x  =  0  (63) 

was  studied  by  Stewartson[31]  in  the  context  of  critical  layers.  The  above  equation  can  be 
solved  using  method  of  characteristics  and  drawing  analogy  with  motion  of  pendulum.  For 
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Figure  8:  Comparison  of  analytical  and  numerical  results  for  the  passive  scalar  problem  at 
X  =  37r/4.  X  is  x(37r/4,  r?,  0)  =  1/(1  +  rj^) 


X  —  7r/4  or  37r/4  one  can  obtain  explicit  expressions  for  x-  At  a;  =  37r/4  if  initially  the 
distribution  of  x  is  x(37r/4,  ry,  0)  =  f{'q)  then  the  solution  at  arbitrary  times  will  be[7], 
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=  / 


/  2  sgn(??) 

V 


dn 


(64) 


where  dn  is  the  Jacobi  elliptic  function  and  K(m)  the  complete  elliptic  integral  of  first 
kind.  Figure  8  shows  the  good  agreement  between  numerical  and  analytical  results  for  a 
Lorentzian  initial  condition.  The  agreement  starts  deteriorating  with  time  due  to  generation 
of  increasing  fine-scales. 


5.3  Holmboe  instabilities 


By  increasing  stratification  we  will  try  to  explore  various  nonlinear  dynamics  for  instability 
scenarios  occurring  in  Holmboe  instabilities  starting  from  stratified  KH  and  proceeding  on 
to  pure  Holmboe  waves.  For  the  purpose  of  exploring  various  nonlinear  states  the  tanh 
model  discussed  before  is  revisited. 
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f  —  2^k  —  1,A  =  1/3,770  =  1  are  the  parameters  maintained  identical  for  all  Holmboe 
simulations.  We  choose  3  values  of  d.  d  =  0.1, 1.5  &  8.  As  visible  from  figure  6,  d  =  0.1 
denotes  stationary  stratified  KH  waves  in  linear  regime,  d  =  1.5,  8  denote  linear  propagating 
Holmboe  waves  (d  =  8  is  not  shown  in  the  figure).  We  will  try  to  find  how  the  nonlinear 
states  behave  for  these  cases.  Does  a  linear  Holmboe  wave  remain  so  always  howsoever 
weak  the  stratification  may  be? 
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5.3.1  Stratified  KH  {d  =  0.1) 

Kelvin-Helmholtz  in  presence  of  weak  stratification  contains  all  the  usual  traits  of  the  un¬ 
stratified  version  -  growth  of  instability  waves,  roll-up  of  the  shear  layer,  vortex  nutation 
(the  horizontal  rocking  of  partially  form  cat’s  eye)  and  finally  culminating  with  the  cele¬ 
brated  cat’s  eye  pattern  or  the  KH  billow.  Besides  the  weakened  growth-rate  due  to  the 
stabilising  effect  of  stratification  the  stratified  version  exhibits  intense  stirring  and  mixing 
of  the  buoyancy  field.  3D  effects  like  elliptical  instability  of  the  billow  and  hyperbolic  in¬ 
stability  of  the  braid  also  emerge  later  in  the  evolution.  Interested  reader  can  find  detailed 
discussions  in  the  article  by  Caulfield  and  Peltier [12]. 

Figure  9  shows  the  evolution  of  buoyancy  (B)  and  vorticity  [Z)  field  at  different  instants 
of  time.  Figure  9(a)-(b)  is  a  snapshot  when  the  flow  exhibits  finite  amplitude  waves  just 
prior  to  the  roll-up.  (c)-(d)  is  at  a  stage  during  the  formation  of  the  billow.  The  flow  has 
departed  from  linear  stability  predictions  and  vortex  nutations  are  observed  at  this  stage. 
Finally  in  figure  9(e)-(f)  the  well  homogenised  KH  billow  is  visible.  Figure  10  shows  the 
time-evolution  of  <  a  measure  of  disturbance  amplification.  Signatures  of  the 

features  discussed  above  can  be  seen  in  the  time-trace. 

5.3.2  Pure  Holmboe  {d  —  8) 

Next  we  explore  a  case  where  the  flow  is  buoyancy  dominant  and  the  shear  layer  dynamics 
are  very  different  from  KH.  For  strongly  stratified  shear  layers  the  flow  exhibits  propagating 
Holmboe  instabilities.  They  are  identified  by  oppositely  travelling  cusps  which  spews  fluid 
from  one  region  to  other.  Compared  to  a  KH  billow  the  Holmboe  waves  are  highly  inefficient 
in  mixing.  The  dynamics  is  very  often  given  by  the  beating  phenomena  of  the  two  Holmboe 
waves.  Laboratory  experiments  on  Holmboe  waves  have  been  done  by  Lawrence  et  al. 
[22],  Pouliquen  et  al.[26]  and  several  researchers  thereafter.  Efforts  towards  using  DNS  to 
understand  Holmboe  wave  dynamics  was  done  by  Smyth  and  co-workers  ([28],  [29]  and 
[30])  amongst  others.  A  recent  work  by  Carpenter  et  al.[9]  considers  both  experimental  and 
numerical  simulation  aspects  of  the  problem. 

From  figure  11  we  observe  that  the  ‘defect’  numerics  are  adept  at  capturing  the  well-known 
traits  of  Holmboe  instability.  Linear  stability  predicts  presence  of  2  unstable  modes  of  same 
growth-rate  but  opposite  phase-speed.  As  expected  the  numerical  solutions  do  exhibit  the 
beating  phenomena  between  these  2  waves.  As  visible  clearly  from  figure  ll(c)-(d)  onwards 
both  the  buoyancy  and  vorticity  field  has  presence  of  elongated  cusp-like  structures  and 
one  can  see  both  waves  before  moving  past  each  other.  Even  at  a  much  later  time  (figure 
ll(e)-(f))  the  cusps  refuse  to  disappear.  If  one  takes  a  look  at  the  evolution  of  disturbance 
amplitude,  <  (figure  12)  the  persistent  beating  pattern  is  immediately  obvious. 

One  cannot  also  help  but  notice  in  figure  11(e)  the  inefficient  nature  of  mixing  of  the 
buoyancy  field  by  Holmboe  waves. 

A  qualitative  comparison  of  the  buoyancy  field  during  Holmboe  instability  computed  using 
‘defect’  numerics  with  those  obtained  from  experiments  and  3D  DNS^[9]  can  be  seen  in 
figure  13.  Though  the  experiments  and  numerics  are  not  in  the  exact  parameter  regime 
as  the  ‘defect’  simulation,  it  is  reaffirming  to  observe  the  present  calculations  ability  to 

^Both  experiments  and  numerics  were  done  for  Reynolds  number,  Re=630,  Richardson  number  being  0.3 
and  Prandtl  number  (temperature  equivalent  to  Sc)  kept  at  700  for  experiment  and  25  for  numerics. 
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(e)  B,t  =  80  (f)  Z,t  =  80 


Figure  9:  Vorticity  and  buoyancy  fields  at  different  instants  of  time  for  stratified  KH  insta¬ 
bility,  /  =  2,  d  =  0.1,  =  1,  A  =  1/3,  r?o  =  1- 
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Figure  10:  Evolution  of  disturbance  amplitude  for  stratified  KH.  q  =  0.9035 
reproduce  the  principal  signature  of  the  shear  layer  with  its  cuspy  arms  as  visible  in  others. 


5.3.3  Holmboe-KH  transition  {d  =  1.5) 

Will  the  persistent  counter-propagating  cusps  be  present  irrespective  of  how  weak  the  strat¬ 
ification  is  as  long  as  linear  stability  predicts  them?  Simulations  for  d—  1.5  as  seen  in  figure 
14  answers  in  negative,  d  =  1.5  is  away  from  the  critical  stratification  {dcrit  —  0.5)  required 
for  transition  from  KH  to  Holmboe  as  predicted  by  linear  stability.  Cusps  get  formed  as 
expected  but  as  visible  in  figure  14(a)-(b)  the  cusps  while  moving  past  each  other  interact 
nonlinearly  and  lock  each  other.  Later  (figure  14(c)- (d))  these  locked  waves  overturn  lead¬ 
ing  to  stirring  of  the  buoyancy  field  leading  towards  a  KH  billow-like  stage  ((e)-(f)).  One 
can  find  reference  to  such  nonlinear  locking  and  subsequent  overt urning/mixing  in  the  ex¬ 
periments  of  Hogg  and  Ivey  [19].  It  is  encouraging  to  see  that  the  reductive  model  obtained 
can  capture  such  intricate  nonlinear  dynamics  at  very  low  computational  cost. 


5.4  Taylor-Caulfield  (TC) 

Instability  due  to  interaction  of  gravity  waves  via  shear  was  one  of  the  intriguing  aspects  of 
Taylor’s  (1931)  analysis.  It  belied  conventional  belief  that  one  needed  an  vorticity  maxima  in 
the  flow  or  unstable  stratification  to  induce  instability  in  stratified  shear  flow.  Moreover  the 
observations  of  staircase  formation  in  shear  flows  provides  a  strong  motivation  to  understand 
how  do  multilayer  density  profiles  interacts  with  shear.  Though  Taylor  introduced  the 
problem  as  a  canonical  stratified  shear  intability  it  wasn’t  until  through  a  series  of  papers 
by  Caulfield  and  co-workers  (linear  stability [10],  experiments [11]  and  DNS[23])  that  a  great 
deal  was  known  about  the  interaction  of  multilayer  density  profiles  with  unbounded  shear 
and  mixing  layer.  The  present  section  will  be  focused  on  unbounded  shear  while  mixing 
layer  will  be  reserved  for  the  next  section. 

The  profile  chosen  for  instability  corresponds  to  the  one  already  dealt  before  in  the  linear 
stability  study  - 


iV2 


d 


sech^ 


(^) 


-h  sech^ 


(^)] 


(67) 


373 


{e)B,t  =  75.5  (f)  Z,t  =  75.5 


Figure  11:  Vorticity  and  buoyancy  fields  at  different  instants  of  time  for  pure  Holmboe 
instability,  f  —  2,d  —  8,k  —  1,  A  —  1/3,  r?o  =  1- 
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Figure  12:  Evolution  of  disturbance  amplitude  for  pure  Holmboe  instability,  q  =  0.4005 


(a)  ‘Defect’  theory 


(b)  Experimental  results 


(c)  3D  DNS 


Figure  13:  Qualitative  comparison  of  buoyancy  field  as  obatined  from  ‘defect’  theory  and 
that  of  both  experiments  and  DNS  (Carpenter  et  al.[9]) 
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(e)  B,t  =  62.5 


X 

(f)  Z,t  =  62.5 


Figure  14:  Vorticity  and  buoyancy  fields  at  different  instants  of  time  for  Holmboe-KH 
transition,  f  —  2,  d  —  1.5,  A:  =  1,  A  =  1/3, 7]o  —  1. 
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Figure  15:  Evolution  of  disturbance  amplitude  for  Holmboe-KH  transition. q  =  0.5537 


for  the  parameters  d  —  5,/c  =  1,A  =  1/8,770  =  2.  The  base-state  has  been  seeded  with 
purely  buoyancy  perturbation  and  vorticity  gets  produced  purely  due  to  baroclinic  effects. 
During  the  linear  instability  stage  the  baroclinically  generated  vorticity  from  the  initial 
conditions  gets  accumulated  at  the  density  interfaces  and  amplifies  (figure  16(a)- (b)).  Once 
the  stationary  TC  waves  get  nonlinearly  saturated  the  buoyancy  forms  a  billow  which  too 
exhibits  nutation  (figure  16(c),  (d)  and  (f)).  The  billow  formation  in  TC  instability  does 
not  involve  the  extensive  overturning  as  observed  in  stratified  KH  but  is  a  more  efficient 
mixing  process  than  pure  Holmboe  waves  (figure  16(e)). 


5.5  TC-Holmboe  transition 


Finally  let  us  consider  a  3-layer  density  profile  embedded  inside  a  mixing  layer.  The  extent 
of  the  density  layer  is  taken  to  be  smaller  than  that  of  the  mixing  layer  (figure  17). 
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Simulation  has  been  carried  out  for  /  =  2,  d  =  4, /c  =  1,  A  =  1/8,771  =  1,772  =  2.  Linear 
stability  predicts  3  unstable  modes  -  a  faster  growing  stationary  TC  mode  and  a  pair  of 
slowly  growing  propagating  Holmboe  wave.  Figure  18(a)- (b)  shows  the  formation  of  a  billow 
in  the  initial  stages  which  starts  spewing  out  cusps  along  its  cusps  ((c)-(d)).  This  cuspy 
billow  persists  and  once  again  a  typical  Holmboe  feature  -  beating  amplitude  pattern  (figure 
19)  -  can  be  seen.  Thus  if  one  were  to  do  an  experiment  on  such  a  profile,  a  propagating 
instability  would  be  observed  despite  linear  stability  making  predictions  for  a  stationary 
billow.  Exploration  of  such  interesting  scenarios  with  great  ease  was  made  possible  by 
virtue  of  working  with  the  ‘defect’  equations. 


6  Nonlinear  stationary  states 

For  stratified  KH  and  TC  instabilities  the  time-trace  for  <  seems  attain  a  sta¬ 

tionary  state  but  for  a  viscous  decay  (figure  10  and  16(f)).  A  qualitative  argument  for  a 
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(e)  Streamwise  averaged  instantaneous  buoy-  (f)  Evolution  of  disturbance  amplitude  for  Taylor- 
ancy  profile  (•  =  J  •  dx)  at  different  instants  Caulfield  instability,  Ci  —  0.5804 
of  time 


Figure  16:  Vorticity  and  buoyancy  fields  at  different  instants  of  time  for  Taylor- Caulfield 
instability,  d  =  5,  fc  =  1,  A  =  1/3, 770  =  2. 
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Figure  17:  Density  layer  embedded  in  a  mixing  layer  (TC-Holmboe  transition) 


stationary  state  also  receives  support  when  one  observes  the  buoyancy  and  vorticity  plots 
at  the  final  time  computationally  attained  (figure  9(f)  and  16(c)).  A  time  evolution  of  the 
plots  show  that  they  are  remarkably  static.  Can  one  provide  quantitative  confirmations  to 
such  guesses? 

Let  us  consider  the  vorticity-buoyancy  defect  equation  under  inviscid,  steady^  assumptions, 


T  —  Bx 

(70) 

V^x  +  ^x^r]  —  0 

(71) 

The  exact  solutions  of  the  above  equations  are. 

B^V{£), 

(72) 

Z=Q{£)+r]V'{£)  =  Q{£)-Br,, 

(73) 

(ti  —  c)^ 

where,  £  = - - - h  ^ 

c  being  the  wave-speed  of  the  state.  V  and  Q  are  arbitrary  functionals  of  the  total  stream- 
function,  £.  Thus  any  nonlinear  stationary  state  needs  to  obey  the  above  functional  re¬ 
lationship.  We  go  ahead  and  try  applying  this  diagnostic  on  both  stratified  KH  and  TC 
instabilities  (c  =  0  in  both  cases). 

At  every  instant  of  time  besides  Z  and  B  the  total  streamfunction,  f,  is  also  computed. 
Thus  at  every  grid  point  (x^,  rjk)  values  oi  B^Z  h£  are  known.  We  construct  pair  of  points 
-  {£ik^  Bik)  and  {£ik^  {Z  -\-  Brj)ik)-  Next  all  the  (f,  B)[{£^  {Z  -\-  Brj))]  pairs  are  plotted  in  the 
B  —  £[{Z  +  Brj)  —  £]  planes.  If  there  indeed  exists  a  functional  relationship  then  the  scatter 
of  points  should  collapse  on  to  a  curve.  That  would  be  a  confirmation  that  the  obtained 
solution  is  indeed  approaching  a  nonlinear  stationary  state. 

First  it  is  applied  on  stratified  KH  (section  5.3.1).  As  visible  from  figure  20  the  initial 
functional  dependence  gets  broken  and  there  is  a  prominent  scatter  during  an  intermediate 
time  -  an  instant  of  strong  instabilities  (figure  20(c)- (d))  -  before  approaching  towards  a 
curve  (figure  20(e)-(f)).  The  near  horizontal  arm(s)  in  both  the  buoyancy  and  vorticity 
plots  corresponds  to  the  uniform  region  outside  the  shear  layer.  From  the  expression  of  the 

the  state  is  moving  with  a  constant  wave-speed  c  then  a  co-moving  reference  frame  is  considered. 
Though  not  presented  here  such  scenario  was  observed  with  simulations  for  asymmetric  stratified  shear 
layer.  One  of  the  Holmboe  waves  is  then  preferred  and  the  system  translates  with  a  constant  speed. 


{e)B,t  =  65.5  (f)  Z,t  =  65.5 


Figure  18:  Vorticity  and  buoyancy  fields  at  different  instants  of  time  for  TC-Holmboe 
transition,  f  —  2,  d  —  A,  k  —  1,  A  —  1/3,  —  l,ri2  —  2. 
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Figure  19:  Evolution  of  disturbance  amplitude  for  TC-Holmboe  transition,  q  =  0.807 


total  stream- function,  one  can  see  why  it  should  be  large  and  negative  outside  the  shear- 
layer.  The  vertical  region,  which  hosts  most  of  the  dynamics  corresponds  to  the  defect. 
Also  the  locations  from  where  points  seem  to  spew  out  vertically,  represent  the  hyperbolic 
points  in  the  flow.  In  the  last  available  buoyancy  plot,  figure  20(e),  the  bulbous  head  is 
representative  of  the  cat’s  eye  region  of  the  flow.  As  visible  from  figure  9(e)  the  B  field  has 
not  homogenised  completely  and  there  is  signature  of  a  dipole-like  structure  in  the  cat’s 
eye,  which  contributes  to  the  bulge. 

As  seen  from  figure  21  similar  observations  are  also  true  for  a  Taylor- Caulfield  case  (section 
5.4). 

Having  identified  the  nonlinear  structures  the  nonlinear  stability  of  the  stationary  states 
can  be  investigated.  If  we  look  for  stationary  points  of  the  energy  under  the  constraints 
that  it  has  to  obey  -  the  conserved  quantities  determined  previously  (section  3)  -  we  have 
the  following  ‘Hamiltonian’  (we  defer  from  calling  the  system  Hamiltonian  since  we  have 
not  checked  for  the  Poisson  bracket  structure  for  it), 

2  /fv  2 

^  V  -rjB-  cqZ  +  F{B)  +  4^  +  ^G{B)  >  (74) 

Zi  z 

where  <  .  >  denotes  average  over  entire  space  as  defined  before.  The  stationary  points  of 
the  system  are, 


SJ{  ^ 

Is 


g  =  G{B),  ^B^  V{S) 
V-F'iB) 


^  Z  — 


G'{B) 


=  Q(^)  +  ^V\E) 


(75) 

(76) 


confirming  what  was  mentioned  in  the  beginning  of  this  section.  Here  5‘K/5Z  is  the  func¬ 
tional  derivative  of  the  functional  fit.  To  check  stability  properties  we  need  to  check  the 
second  variation  fit. 


=  <  -5^5Z  +  +  2G'{B)SZSB  > 

=  <{SZ  +  g'(B)SB}^  +  -  (SZfj  +  {F"(B)  +  Zg"(B)  -  (g'(B)f}(SBf  > 
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(c)  t  =  14.5 


(d)  t  =  14.5 


Figure  20:  Nonlinear  stationary  state  diagnostics  for  stratified  KH  instability.  Time  evolu¬ 
tion  of  the  dependence  of  B  with  £  is  shown  in  (a),  (c)  and  (e)  and  that  of  Z  —  rjB'  with  £ 
in  (b),  (d)  and  (f). 


In  this  form  a  comment  cannot  be  made  about  nonlinear  stability  as  5‘^'K  is  indefinite.  But 
if  J  dxdri{—5ip5Z  —  >  0  is  positive  definite  then  the  condition  for  stability  will  be  - 

+  Zg"{B)  -  {G\B)f}  >  0 


Now  the  question  is  under  what  cases  is  the  assumption,  J  dxdri{—6^p6Z  —  (dZ)^^}  >  0,  true. 
In  ‘defect’  theory,  for  an  unbounded  domain,  ^  and  Z  can  be  connected  in  Fourier  space 
as. 


~  j  Z{k,r],t)dr] 


(77) 


Thus  using  Parseval’s  theorem  we  have 


J  dxdrji-S^SZ  -  (5Z)’^}  >  0 

J  dkd7]{-SipSZ*  -  \SZ\‘^}  >  0 


r  1 

f 

) 

\  2k 

/  drjSZ 

-  /  drjlSZl"^  }  >  0 


(78) 

(79) 

(80) 


Presently  a  detailed  study  in  the  above  aspects  is  lacking  but  it  is  hypothesized  that  there 
would  exist  a  band  of  wavenumbers  in  which  a  stratified  shear  flow,  modelled  by  ‘defect’ 
theory,  will  be  stable.  For  higher  wavenumbers  one  can  always  find  instability. 


7  Conclusion 

Exploiting  the  sharp  variations  present  in  stratified  shear  layers  which  occurs  over  small 
regions,  ‘defect’  vorticity-buoyancy  equations  have  been  derived.  This  equation  maintains 
all  the  necessary  physics  -  unsteadiness,  nonlinearity  and  viscosity.  Linear  stability  cal¬ 
culations  using  the  ‘defect’  equations  becomes  significantly  convenient  with  the  present 
equations.  But  the  big  strength  of  ‘defect’  equations  was  found  to  be  in  doing  nonlinear 
calculations  for  various  stratified  shear  layer  profiles.  Unlike  conventional  numerics  done 
for  Bousinessq  Navier-Stokes  equations  the  ones  for  ‘defect’  equations  can  be  done  with 
minimum  computational  effort. 

We  have  been  able  to  explore  some  interesting  cases  concerning  nonlinear  evolution  of  strat¬ 
ified  shear  layers.  Though  a  great  deal  is  known  about  the  stratified  KH  instability,  we  were 
able  to  address  some  intriguing  questions  regarding  the  transition  from  KH  instability  to 
Holmboe  waves  and  also  whether  a  Holmboe  instability  persists  nonlinearly  with  decreasing 
stratification.  Very  little  is  known  about  the  nonlinear  dynamics  of  Taylor- Caulfield  insta¬ 
bilities.  An  effort  has  been  made  to  probe  TC  instability  in  greater  detail  while  working 
with  ‘defect’  equations.  Coexistence  of  a  TC  billow  with  a  Holmboe  cusp  was  also  observed 
while  studying  density  layers  embedded  in  mixing  layer. 

A  diagnostics  was  used  with  success  to  identify  nonlinear  stationary  structures  in  stratified 
KH  and  TC  instabilities.  Motivated  by  visual  arguments  a  check  was  done  on  the  buoyancy, 
vorticity  and  total  streamfunction  field  to  determine  if  the  flow  is  attaining  an  invariant 
form.  Efforts  to  perform  nonlinear  stability  have  so  far  been  inconclusive. 
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Figure  22:  Time  evolution  of  the  scalar  variance 


Being  built  for  flow  with  sharp  interfaces  or  local  distortions,  the  ‘defect’  theory  has  immense 
potential  for  stratified  shear  layer  studies.  TC  instabilities  is  far  from  being  well-understood 
and  it  will  be  very  interesting  to  explore  a  wide  range  of  parameter  space  armed  with  ‘defect’ 
theory.  An  aspect  completely  overlooked  in  the  present  analysis  is  the  estimation  of  mixing 
in  the  considered  flows.  Since  Sc  ^  oo  was  the  domain  of  interest  an  estimate  of  stirring 
was  made  based  on  the  scalar  variance,  <  iVBp  >.  A  low  value  of  <  >  denote 

homegenisation  of  the  scalar  field.  Figure  22  does  indeed  highlight  an  already  discussed 
aspect,  based  on  visual  grounds,  of  mixing/stirring  efficiency  being  maximum  for  stratified 
KH,  followed  by  TC  and  then  Holmboe.  A  more  careful  analysis  needs  to  be  done  in  this 
regard  and  comparisons  made  with  existing  DNS  and  experimental  results. 

The  ‘defect’  equations  are  amenable  to  analytical  approaches  more  easily  than  the  govern¬ 
ing  equations  from  which  they  have  been  derived.  We  haven’t  really  used  this  to  obtain 
the  solutions  in  the  nonlinear  states  in  presence  of  dissipation.  A  quasi-steady  critical  layer 
theory  approach  is  being  investigated  under  the  framework  of  ‘defect’  equations  to  answer 
the  asymptotic  form  of  the  nonlinear  stationary  states.  Finally  a  potential  future  work 
which  has  been  left  out  from  the  present  study  is  the  study  on  the  viscous  spectrum,  issue 
of  quasimodes  and  transient  growth. 
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